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Abstract 


During  the  period  of  sponsorship  under  the  JSEP  fellowship,  the  researcher  has  studied  tomographic 
imaging,  the  mathematical  process  that  underlies  a  number  of  technologies  such  as  Synthetic  Aper¬ 
ture  Radar,  X-ray  computed  tomography,  Magnetic  Resonance  Imaging,  etc.  In  particular,  his 
research  has  focused  on  two  key  aspects  of  tomography.  The  first  is  the  study  of  tomography  in 
the  presence  of  motion.  In  this  area,  he  was  the  first  to  establish  conditions  for  unique  and  stable 
reconstruction  in  the  presence  of  rigid  body  motion.  The  second  aspect  is  that  of  fast  algorithms 
for  tomography.  He  has  developed  new  advanced  algorithms  that  permit  tomographic  imaging  to 
be  carried  out  much  faster  than  was  previously  possible.  In  the  process  of  studying  these  two  key 
areas,  he  has  developed  new  and  fundamental  results  in  signal  processing. 
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Chapter  1 

Introduction 

Goals: 

Our  research  work  focused  on  two  goals.  Our  first  goal  was  to  study  tomographic  reconstruction 
of  time- varying  densities  by  means  of  motion  models.  Our  second  goal  was  to  study  and  develop  a 
new  class  of  fast  algorithms  for  tomography. 

Motivation: 


Tomography  is  a  technique  by  which  a  function  defined  in  the  plane  is  reconstructed  from 
a  collection  of  line  integral  measurements  at  different  orientations  [1],  It  has  found  numerous 
applications  in  the  fields  of  medical  imaging,  nondestructive  evaluation,  radio  astronomy,  Synthetic 
Aperture  Radar  (SAR),  electron  microscopy  based  tomography,  and  nondestructive  evaluation 
(NDE)  among  others.  As  such,  it  is  a  pervasive  technique  and  has  been  well  studied. 

The  mathematical  characterization  of  the  tomographic  problem  was  initially  described  by  Radon 
in  1917  [2],  Radon  constructed  an  inversion  formula  for  tomographic  data  in  2D,  when  the  line 
integrals  were  available  from  all  possible  directions.  With  the  introduction  of  the  Computed  To¬ 
mography  (CT)  scanner,  it  became  possible  to  collect  these  line  integral  measurements,  albeit  at 
a  finite  set  of  directions.  Much  work  over  the  last  few  decades  has  focused  on  reconstruction  tech¬ 
niques  from  such  finite  sets  of  data  [3].  Many  techniques  have  been  described  and  implemented, 
and  tomographic  reconstruction  is  generally  considered  a  solved  problem. 
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However,  as  imaging  technology  improves,  the  limitations  of  the  traditional  tomographic  for¬ 
mulation  are  beginning  to  manifest  themselves  [4],  We  have  studied  two  of  these  limitations  in  this 
work.  The  first  is  the  limitation  of  assuming  a  temporally-static  2-D  density.  In  most  practical 
situations,  the  acquisition  is  time-sequential,  so  that  different  sets  of  measurements  are  taken  at 
sequential  moments  in  time.  If  the  2-D  density  being  imaged  is  evolving  in  time,  then  the  measure¬ 
ments  violate  the  assumption  that  each  measurement  is  of  the  same  2-D  density  This  violation 
leads  to  imaging  artifacts  that  can  render  the  reconstructions  useless.  Indeed,  extensive  documen¬ 
tation  exists  of  the  types  of  artifacts  that  occur  in  medical  imaging  scenarios  as  a  result  of  voluntary 
(movement  of  the  arms,  head,  etc.)  and  involuntary  patient  motion  (cardiac,  gastrointestinal,  etc), 
and  the  resulting  effect  on  the  diagnostic  utility  of  the  reconstruction  [5, 6].  The  second  limitation 
we  have  examined  is  that  of  computational  complexity.  Traditional  techniques  for  reconstruction 
are  generally  computationally  demanding.  As  real-time  reconstruction  becomes  increasingly  impor¬ 
tant.  the  push  for  faster  reconstruction  algorithms  with  little  distortion  becomes  more  significant 
as  well  [4].  With  the  recent  ground-breaking  introduction  of  a  hierarchical  reprojection  algorithm 
in  [7],  it  may  be  possible  to  achieve  this  goal. 

Organization: 

In  Chapter  2,  we  review  the  state  of  the  art  on  the  issues  of  tomography  in  the  presence  of 
motion  and  on  fast  algorithms  for  tomography.  In  Chapter  3,  we  describe  the  results  achieved 
while  sponsored  by  the  JSEP  fellowship.  This  includes  fundamental,  new  results  in  tomographic 
reconstruction  in  the  presence  of  motion,  as  well  as  a  new  class  of  fast  algorithms  for  tomography. 
Some  fundamental  results  in  signal  processing  that  we  developed  as  a  consequence  of  our  primary 
research  interests  are  also  described  there.  In  Chapter  4,  we  discuss  our  conclusions  regarding 
the  work  we  have  completed.  We  have  attached  a  list  of  submitted  and  published  papers  written 
during  the  duration  of  this  fellowship.  In  accordance  with  the  requirements  for  this  final  report, 
those  papers  are  not  attached  as  appendices. 
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Chapter  2 

Review 

i  I.  Motion  in  Tomography 

A.  Formulation 

As  traditionally  formulated,  the  tomographic  reconstruction  problem  is  to  reconstruct  an  image 
supported  on  the  unit  disk  82,  from  a  set  of  P  parallel  beam  projections  at  view  angles  6t  G  [0, 7r). 
If  we  assume  that  the  image  belongs  to  the  Hilbert  space  of  square-integrable  functions  supported 
01182  (denoted  £2(82)),  then  we  can  define  a  Radon  transform  operator  1Z  :  £2(82)  l~>  ^2*(®i)  85 

where  1^(8 1)  ls  the  P-wise  Cartesian  product  of  £2(81),  the  Hilbert  space  of  square  integrable 
function  supported  on  Bi  =  [—1,1].  The  traditional  motion-free  tomographic  reconstruction  prob¬ 
lem  is  to: 

Definition  1  (Motion-Free  Tomography)  Estimate  f  from  noisy  measurements  oflZf. 

Although  this  problem  is  well-studied,  its  formulation  ignores  one  important  physical  consid¬ 
eration.  In  most  data  acquisition  scenarios,  these  measurements  are  obtained  sequentially ,  and  / 
may  change  with  time.  When  this  change  is  significant,  the  motion-free  tomographic  model  is  no 
longer  valid,  and  artifacts  can  occur  in  the  reconstruction. 

To  handle  the  dynamic  nature  of  /,  we  must  redefine  the  problem.  We  now  consider  /  to  be 
a  square  integrable  3D  density,  where  the  third  dimension  is  a  time  dimension.  It  is  supported  on 


TZf(s,i)  =-  /  f(s  cos  6i  +  r  sin  Qi,  s  sin  6t  —  r  cos  #*)  dr, 
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the  cylinder  F3  =  82  x  [0, 1],  and  is  an  element  of  the  Hilbert  space  of  such  functions,  denoted 
L2(¥s).  To  convert  the  motion-free  formulation  to  a  more  general  one,  we  consider  the  following 

‘  i  . 

time-sampled  Radon  Transform :  1Zt  :  Z2OF3)  Bi) 

i)  =  J  f(scos0i  +r  sin  6^$  sin  6i  —  r  cos  0^7$)  dr, 

where  T{  G  [0,1].  Note  that  this  formulation  allows  for  multiple  projections  to  be  acquired  at 
the  same  time  instant  Tj,  simply  by  choosing  =  tj  for  each  pair  of  view  angles  i,j  acquired 
simultaneously.  The  vector  6  of  view  angles  and  the  vector  r  of  time  instants  jointly  characterize  the 
imaging  geometry.  If,  for  example,  the  imaging  device  is  an  axial  CT  scanner,  then  the  acquisition 
geometry  is  uniform ,  meaning  that  0*.  =  A;A0,  t*.  =  kT .  This  acquisition  geometry  also  occurs  in 
the  tomographic  formulation  of  SAR. 

The  problem  to  be  solved  now  is  the  following  one: 

Definition  2  (Time-sampled  Tomography)  Estimate  f  from  noisy  measurements  of  7Ztf  • 

B.  Past  Research 

A  diverse  set  of  results  have  already  been  obtained  on  the  time-sampled  tomography  problem. 
Broadly,  these  results  can  be  divided  into  two  parts.  The  first  deal  primarily  with  the  parameters 
6  and  r,  and  how  they  can  be  modified  to  make  solution  of  the  time-sampled  problem  feasible. 
The  second  set  of  results  assume  a  fixed  choice  of  6  and  r  (usually  a  uniform  geometry),  and  then 
study  how  different  motion  models  affect  solvability  of  the  time-sampled  tomography  problem. 

Modified  Acquisition 

Several  researchers  have  looked  at  different  means  of  choosing  the  geometry  so  as  to  best  solve 
the  time-sampled  tomography  problem.  The  definitive  work  in  this  area  is  to  be  found  in  [8,9]. 
There,  the  authors  examine  the  case  in  which  /  is  spatially  and  temporally  bandlimited.  They 
construct,  using  an  elegant  lattice-packing  argument,  the  optimal  choices  for  6  and  r  under  these 
circumstances,  as  well  as  bounds  on  the  minimum  angular  and  temporal  rates  necessary  to  avoid 
aliasing  in  the  reconstruction. 

Similar  results  were  obtained  in  [10],  in  which  methods  are  described  which  adaptively  choose 
^optimal"  6  based  on  /  itself.  There,  the  authors  also  provide  efficient  reconstruction  methods 
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using  the  wavelet  localization  of  the  Radon  transform  [11].  The  results  provided  are  encourag¬ 
ing,  demonstrating  that  the  time-sampled  tomography  problem  is  at  least  theoretically  solvable, 
provided  the  view  angles  can  be  placed  correctly. 

An  entirely  different  approach  in  medical  imaging  is  that  of  so-called  gated  acquisition .  In  this 
approach,  it  is  generally  assumed  that  /  is  periodic,  and  that  some  indicator  of  the  phase  within 
the  period  can  be  measured1.  Specific  approaches  to  gating  are  described  in  [12-14].  Generally, 
periodicity  is  exploited  to  reduce  the  time  varying  problem  into  a  collection  of  static  problems. 
Gating  approaches  suffer  from  both  a  requirement  of  cooperation  on  the  part  of  the  patient,  and 
an  assumption  of  “proper”  biological  operation  to  guarantee  periodicity.  Thus,  imaging  under 
unusual  medical  circumstances  (such  as  cardiac  arrythmia)  can  defeat  gating  based  approaches. 
Furthermore,  gated  approaches  are  clearly  useless  for  imaging  nonperiodic  phenomena.  Tekalp 
offered  artifact  correction  methods  for  periodic  models  without  gating,  in  which  the  periodicity 
was  extracted  from  the  MRI  data  [15].  There,  the  out-of-slice  motion  was  modeled  as  a  modulation 
of  the  projection  data.  Haacke  et  al.  have  proposed  Projection  Onto  Convex  Sets  (POCS)  based 
reconstruction  methods  for  periodic  motion  models,  when  the  motion  is  partially  known  [16]. 

The  general  difficulty  of  these  approaches  in  general,  is  that  they  require  some  modification  of 
the  method  by  which  data  is  acquired.  In  many  cases  of  practical  interest,  such  as  SAR,  NDE, 
and  Spiral  CT,  this  is  impractical  or  impossible.  Although  some  technologies  such  as  MRI  or 
the  Electron  Beam  Tomography  [17]  scanner  permit  great  flexibility  in  the  acquisition  geometry, 
in  the  majority  of  applications,  the  geometry  is  fixed  by  physical  considerations.  For  example, 
in  spotlight-mode  SAR  the  view  angles  are  determined  by  the  relative  positions  of  the  antenna 
and  the  target  patch.  Thus,  physical  constraints  on  the  motion  of  the  antenna  (which  might,  for 
example,  be  attached  to  an  aircraft)  limit  the  flexibility  of  the  acquisition  sequence.  Likewise,  in  CT 
scanners,  the  radiation  source  is  typically  mounted  on  a  high  speed  rotating  ring.  The  rotational 
motion  restricts  the  acquisition  geometry  to  be  uniform  2. 

Model  based  approaches 

Other  results  on  time-sampled  tomographic  imaging  generally  fix  the  imaging  geometry,  and 

Examples  of  indicators  include  ECG  signals  for  cardiac  imaging  and  chest  displacement  for  imaging  of  the  lungs. 

2 It  is  possible  in  theory  to  change  the  angular  sampling  pattern  of  a  CT  scanner  by  turning  the  x-ray  source  on 
and  off  at  the  appropriate  times. 
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then  study  how  different  models  for  the  time  evolution  of  /  affect  solvability.  Different  researchers 
have  adopted  different  models,  generally  motivated  by  physical  concerns.  For  respiratory  compen¬ 
sation  in  tomographic  imaging,  global  and  local  affine  models  have  been  studied  [18, 19].  For  the 
global  affine  model,  the  function  f(x,t)  is  related  to  some  unknown  static  function  fo{x)  by 

f{x,Ti)  =  f0(DiX  +  bt,t) 


where  D{  E  R2x2  are  known  diagonal  matrices,  and  b*  E  R2  are  known  shifts3.  The  local  affine 
model  is  similar,  except  that  the  affine  model  is  applied  locally. 

Some  simple  translational  motion  models  have  been  studied  in  the  context  of  different  geome¬ 
tries,  including  2D  MRI4  [20-22],  SAR  [23],  and  helical-CT  [24].  These  model-based  approaches 
have  been  fairly  successful,  provided  that  the  models  adequately  describe  the  motion. 

Another  class  of  models  have  been  proposed  by  Liang  et  al  [25-27].  Here,  a  'priori  knowledge 
of  the  spatial  distribution  of  /  for  each  sample  instant  is  incorporated  in  the  form  of  a  specially 
chosen  set  of  basis  functions.  A  separable  model  of  the  form 


/(x,  i)  —  ^  ^  Cn(^)<ftn fo) 

n 


is  then  used  to  model  the  temporal  evolution  of  /.  Assuming  that  the  basis  functions  4>n(x) 
can  be  chosen  appropriately,  these  models  can  permit  the  study  of  ‘'compartmentalized”  dynamic 
changes,  such  as  may  occur  in  spectroscopy  or  in  contrast  studies.  The  coefficients  Cn(t )  are  uniquely 
determined  by  the  data,  provided  that  there  are  a  sufficiently  small  number  of  degrees  of  freedom 
in  the  model  [25]. 

In  general,  motion-model-based  approaches  permit  the  incorporation  of  specific  prior  informa¬ 
tion  into  the  behavior  of  /.  This  can,  of  course,  also  be  their  Achilles  heal,  particularly  if  the 
resulting  reconstruction  technique  is  highly  sensitive  to  deviations  from  the  model. 


*The  diagonal  matrices  and  shifts  were  obtained  through  external  measurements,  similar  to  gating,  although 
periodicity  was  not  assumed. 

4 Although  the  cited  works  were  done  using  the  rectilinear  sampling  strategy,  and  not  the  projection  mode  of 
acquisition,  which  is  tomographic. 


II.  Fast  Algorithms 

A.  Problem  Formulation 

In  the  realm  of  motion-free  tomography,  there  are  difficulties  that  arise  in  practice  in  the 
routine  solution  of  the  tomographic  reconstruction  problem.  One  important  problem  is  the  high 
computational  requirement  that  real-time  tomographic  reconstruction  places  on  existing  computers. 
As  medical  X-ray  CT  scanners  become  more  sophisticated,  the  reconstruction  process  becomes  a 
major  bottleneck  in  the  processing  of  data  for  visualization  [4].  Hence,  there  is  a  growing  need 
for  fast  reconstruction  algorithms  that  can  produce  diagnostic  quality  reconstructions  in  a  small 
fraction  of  the  amount  of  time  currently  needed. 

B.  Review 

The  traditional  approach  to  accelerating  the  reconstruction  process  is  to  employ  special  purpose 
hardware,  in  the  form  of  pipelined  and  hardwired  architectures  that  perform  the  backprojection 
operation  at  high  speed  [28-30].  Parallel  computing  is  another  method  used  to  speed  up  backpro- 
jection.  These  traditional  approaches  suffer  from  two  major  problems.  First,  the  special  purpose 
hardware  quickly  grows  obsolete,  as  general  purpose  microprocessors  become  increasingly  faster. 
Second,  they  only  offer  0(1)  (i.e.  constant  factor)  speedups  over  a  standard  implementation  on  a 
serial  processor. 

Recently,  a  number  of  fast  reconstruction  algorithms  have  been  proposed  [31-33].  Each  of  these 
fast  algorithms  offers  0(N/  log  N)  speedup  over  direct  implementation,  where  N  is  the  linear  size 
of  the  reconstructed  image  in  pixels.  The  advantages  of  such  algorithms  over  direct  implementation 
for  large  N  are  obvious.  With  current  sub-millimeter  resolution  CT  scanners,  N  =  103  and  larger 
images  are  becoming  common  place.  Fast  algorithms  promise  one  to  two  orders  of  magnitude 
speedup  for  such  images.  And  yet,  the  use  of  fast  algorithms  is  not  common.  The  main  reason 
appears  to  be  that  each  of  these  fast  algorithms  suffers  from  some  fundamental  limitation  that 
leads  to  overall  poor  performance. 

The  oldest  fast  algorithms  (and  perhaps  the  most  obvious)  are  known  as  Fourier  Reconstruction 
Algorithms  [34-38]  (FRA).  Based  on  the  Fourier  Slice  Theorem,  these  algorithms  interpolate  the 
tomographic  data  from  the  polar  on  which  they  are  acquired  to  a  rectangular  grid,  both  in  the 
Fourier  domain  (see  Figure  2.1.  An  inverse  2D-FFT  can  then  be  used  to  recover  the  cross-sectional 
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Figure  2.1:  FRAs  interpolate  the  Fourier  data  from  the  polar  grid  of  the  tomographic  measurements, 
to  the  Cartesian,  FFT  grid. 

density.  In  theory,  this  yields  an  0(N2  log  N)  algorithm  for  reconstruction.  But  the  reports  of 
the  performance  of  the  FRA  have  been  disappointing  at  best  [39,40].  In  general,  the  precision 
necessary  in  the  interpolation  step  results  in  an  algorithm  that,  for  reasonable  sized  images,  is  only 
slightly  faster  than  the  straightforward  reconstruction  algorithm,  which  requires  0(N3)  operations 
[41].  The  latest  incarnation  of  the  FRA  is  known  as  gridding  [31,42],  and  includes  all  other  FRA’s 
*ts  special  cases.  Gridding  provides  good  reconstructions,  but  only  provides  small  speedups  for 
reasonable  image  sizes. 

Next  in  line  are  the  linogram  reconstruction  algorithms  [32].  These  algorithms  are  based  on 
the  assumption  that  the  view  angles  lie  at  certain  special  angles,  so  that  when  the  projections  are 
cleverly  resampled,  the  reconstruction  can  be. performed  using  a  series  of  ID  inverse  FFTs.  The 
requirements  on  the  angles  (and  the  radial  resampling  of  the  projections  necessary  in  practice) 
make  the  resulting  algorithms  unattractive  from  a  practical  standpoint. 

The  third  type  of  fast  reconstruction  algorithm  is  a  multiscale  approach  [33].  In  this  algorithm, 
a  sequence  of  nonuniform  sample  “grids”  are  used  to  merge  the  projection  data  in  progressive 
steps,  until  the  samples  represent  the  reconstruction  itself.  Unfortunately,  this  algorithm  generates 
blurred  reconstructions,  which  must  then  be  corrected  using  deconvolution  techniques.  As  this 
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algorithm  is  fairly  new,  less  is  available  as  to  its  performance  characteristics. 

In  addition  to  fast  algorithms  for  reconstruction,  there  is  also  a  need  for  fast  algorithms  to 
perform  a  variety  of  tomography-related  computations.  All  of  the  methods  described  so  far  can  be 
used  for  fast  backprojection.  Also  of  use  are  fast  reprojection  methods,  in  which  projection-data  are 
synthesized  from  an  image.  When  fast  reprojection  and  backprojection  algorithms  are  combined, 
they  can  form  the  basis  of  fast  and  powerful  iterative  reconstruction  techniques.  In  [43],  just 
such  an  iterative  technique  is  proposed.  There,  the  reprojection  and  backprojection  operations  are 
combined  into  a  single  operation  that  can  be  computed  efficiently  using  2D  FFT’s.  The  resulting 
technique  is  both  fast  and  powerful.  Unfortunately,  only  operations  on  the  measurement  data 
that  are  radially  shift  invariant  can  be  inserted  between  the  reprojection  and  backprojection  steps 
without  destroying  the  convolutional  nature  of  the  algorithm.  This  limitation  is  rather  restrictive, 
and  significantly  reduces  the  usefulness  of  the  results  in  [43]. 

Recently,  a  fast  reprojection  algorithm  was  proposed  that  appears  to  suffer  few  of  the  problems 
that  plague  the  other  fast  algorithms  [7].  This  fast  algorithm  is  based  on  a  fundamentally  new 
decomposition  of  the  Radon  transform.  It  is  capable  of  computing  0(N )  projections  from  an  N  x  N 
image  in  0(N2  log2  N)  time,  with  small  distortion.  Our  work  on  fast  algorithms  in  tomography  are 
based  on  the  decomposition  described  in  [7]. 
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Chapter  3 

Work  Completed 

I.  Rigid  Body  Motion 

Our  work  on  tomography  and  motion  has  focused  primarily  on  the  use  of  specific  kinds  of 
parametric  motion  models  to  solve  the  time-sampled  tomography  problem  with  a  fixed  acquisition 
geometry.  A  surprising  result  that  was  first  proposed  simultaneously  by  Goncharov  and  Salzman, 
is  that  if  /  is  defined  by 

f(x,Tt)  =  fo(QiX  +  bi) 

where  Q{  is  a  unitary  matrix  in  R2x2,  and  6;  is  sufficiently  small  so  that  f(x,  r i)  is  supported  on  ®2, 
then  the  parameters  and  6*  that  describe  the  “motion”  at  each  sample  instant  can  be  determined 
from  the  projection  data  alone  [44, 45].  Both  Goncharov  and  Salzman  were  primarily  interested  in 
these  results  for  the  purposes  of  solving  the  alignment  problem  in  Cryo-Electron  Microscopy.  This 
problem,  illustrated  in  Figure  3.1,  is  to  determine  the  orientations  of  a  number  of  identical  particles 
embedded  in  a  solid,  from  a  single  projection,  or  very  small  number  of  projections.  As  discussed  in 
[45].  the  problems  are  substantially  different  in  2D  and  3D  (as  are  many  aspects  of  tomography), 
and  the  problem  is  nearly  trivial  in  3D. 

A.  Uniqueness 

Although  the  idea  that  tomography  in  the  presence  of  rigid  body  motion  was  suggested  by 
Salzman  and  Goncharov,  they  failed  to  demonstrate  rigorously  that  the  parameters  of  the  rigid 
body  motion  model  were  indeed  uniquely  recoverable.  A  guarantee  of  uniqueness  may  be  critical  in 
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Figure  3.1:  How  the  unknown  view  angle  problem  arises  in  Cryo-EM  imaging  of  macromolecules. 

sensitive  applications,  such  as  medical  imaging,  where  nonunique  solution  may  lead  to  misdiagnosis. 

A  major  accomplishment  of  our  research  under  this  fellowship  was  to  demonstrate  carefully 
and  rigorously  that  the  motion  parameters  are  uniquely  determined  by  the  projection  data,  and  to 
exactly  determine  the  conditions  that  must  be  satisfied  for  these  uniqueness  conditions  to  hold.  In 
[46.47],  we  rigorously  demonstrated  the  following  fact: 

Result  1  For  almost  all  /o,  and  almost  all  6  with  P  >  24,  Qx  and  b{  are  uniquely  determined  by 

ntf. 

Furthermore,  our  techniques  laid  the  foundation  for  future  work,  in  the  form  of  more  complex  mo¬ 
tion  models,  and  a  mechanism  for  determining  conditions  for  unique  recovery  of  motion  parameters 
from  projection  data. 

B.  Stability 

Having  established  that  the  motion  model  parameters  were  uniquely  determined  from  the  pro¬ 
jection  data,  the  next  step  was  to  study  the  effect  of  noise  in  the  projection  data  on  the  solution 


for  the  motion  parameters.  In  [48],  we  proved  the  following  fact: 

Result  2  For  almost  all  /o,  and  almost  all  6  with  P  >  24;  Qx  and  bi  are  uniquely  determined  by 
IZtf,  and  furthermore,  the  problem  of  determining  and  bi  from  IZtf  in  the  presence  of  noise  is 
stochastically  stable,  in  the  sense  that  the  Cramer-Rao  bounds  are  finite . 

We  also  elucidated  the  connection  between  the  Cramer-Rao  bounds  and  a  deterministic  notion 

of  stability  in  [49].  With  the  aid  of  this  analysis,  we  obtained  the  following  additional  result: 

'  •  .  '  .  .  ' 

Result  3  For  almost  all  fo,  and  almost  all  6  with  P  >  24,  Qx  and  bi  are  uniquely  determined  by 
7ltf,  and  furthermore ,  the  problem  of  determining  Qx  and  b{  from  TZtf  in  the  presence  of  noise  is 
deterministically  stable,  in  the  sense  that  the  least  squares  estimate  depends  smoothly  on  the  data 
for  any  sufficiently  small  perturbation  of  the  data. 

Showing  that  the  bounds  were  finite  was  a  first  step.  We  also  obtained  the  following  (numerical) 
result: 

Result  4  The  rotational  and  translational  components  of  the  rigid  motion  can  (according  to  the 
Cramer-Rao  bounds)  be  estimated  to  within  a  fraction  of  a  degree  and  a  small  percentage,  respec¬ 
tively,  at  an  SNR  of  30dB. 

These  results  suggest  that  the  problem  of  determining  Qx  and  bt  from  !Ztf  is  feasible,  provided  the 
noise  level  is  sufficiently  small. 

Our  final  result  on  this  problem  (also  reported  in  [48])  was  a  simple  algorithm  to  estimate  Q{ 
and  bi'. 

Result  5  A  simple  and  efficient  algorithm  exists  that  can  estimate  Qx  and  bl7  and  that  achieves 
the  Cramer-Rao  bound  for  SNR  >  20dD  on  a  computer  phantom. 

It  is  important  to  point  out  that  while  none  of  the  previous  results  made  assumptions  as  to  the 
nature  of  Qz  and  the  algorithm  we  proposed  made  some  assumptions.  Specifically,  the  proposed 
algorithm  assumes  that  there  exists  a  permutation  7 r  of  the  projections  such  that  ~  Q7r(i+1) 

and  b^i)  ~  Thus,  under  some  (unknown)  reordering  of  the  samples,  the  motion  described 

by  Q  and  b  was  small  between  samples.  The  algorithm  then  searches  for  the  correct  permutation 
of  the  projections.  In  [48],  we  demonstrated  that  this  strategy  works  well  if  a  large  number  of 
projections  are  available  without  significant  “gaps”  between  them. 
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II.  Fast  Algorithms  for  Tomography 

A.  Fast  Backprojection 

As  we  mentioned  in  Chapter  2,  the  introduction  in  [7]  of  a  fundamentally  new  decomposition  of 
the  Radon  transform  has  lead  to  a  new  class  of  fast  algorithms  for  reprojection,  i.e.  the  computation 
of  projections  from  a  discrete  image.  Unlike  direct  computational  methods,  which  require  0(N 3) 
operations  to  compute  N  projections  from  and  N  xN  image,  the  new  fast  algorithm  of  [7j  promises 
the  same  results,  but  in  0(N 2  log2  N )  operations.  For  large  image  sizes,  the  potential  image  speedup 
is  enormous. 

The  operation  of  reprojection  is  useful  in  tomographic  reconstruction,  particularly  in  the  context 
of  iterative  algorithms.  But  the  operation  of  backprojection  is  even  more  prevalent.  Given  a 
collection  of  P  projections  gc(r,  p )  for  r  E  R.  and  p  G  {1, . . .  ,  P),  backprojection  requires  computing 

P 

f(hj)  =  9c(i  cos  6p  +  j  sin  ep,P)  (2.i) 

P= 1 

where  6p  is  the  p-th  view  angle.  A  straightforward  discretization  of  (2.1)  leads  to  an  0(N 3) 
algorithm  if  P  =  O(N). 

However,  we  have  been  able  to  use  the  same  type  of  decomposition  described  in  [7]  to  generate 
a  fast  algorithm  to  perform  this  same  computation  in  0(jV2  log 2  N)  operations  [50]: 

Result  6  A  simple >  fast  and  controllable  algorithm  to  compute  the  backprojection  of  P  projections 
oil  an  X  x  N  pixel  grid  in  0(NP  log2  iV)  time. 

Experimental  results  demonstrate  that  the  new  algorithm  achieves  the  promised  iV/log2  N  speedup 
for  large  images: 

Result  7  The  resulting  fast  algorithm  can  compute  backprojections  for  1024  x  1024  pixel  images 
ov(  r  00  times  faster  than  direct  methods ,  and  faster,  and  with  lower  distortion  than  any  other  fast 
rt  construction  method  available  at  the  time. 

III.  Fundamental  Results  in  Signal  Processing 

As  a  byproduct  of  our  primary  research  goals,  we  have  developed  some  fundamental  results  in 
signal  processing. 
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A.  Global  Bounds  on  Estimation 

In  many  signal  and  image  processing  problems,  one  must  estimate  a  vector  parameter  0  given 
a  noisy  message  g(t ,  0)  that  is  parameterized  by  6.  Furthermore,  in  many  situations  of  practical 
interest,  the  g  depends  on  some  components  of  6  in  a  periodic  manner.  Examples  of  problems 
of  this  kind  include  phase  estimation  from  a  periodic  pulse  train.  In  this  example,  the  message 
g  depends  on  the  phase  in  a  periodic  manner.  Image  registration  is  another  example.  Here,  the 
vector  6  may  represent  the  rotation  and  shift  necessary  to  co-align  two  images.  In  that  case,  the 
dependency  of  g  on  the  rotational  parameter  is  periodic.  Another  example  is  orientation  estimation. 
In  this  example,  the  g{t ,  6)  represents  a  signal,  like  a  radar  or  sonar  return  of  an  object  in  unknown 
position.  The  dependence  of  g  on  the  positional  parameters  corresponding  to  roll,  pitch  and  yaw, 
is  periodic. 

In  all  of  these  cases,  some  of  the  parameters  affect  the  message  in  a  periodic  manner,  and  the 
other  parameters  affect  the  message  in  a  nonperiodic  manner.  For  parameters  of  the  first  kind,  the 
correct  choice  for  measuring  distortion  is  to  use  a  periodic  distortion  measure.  For  example  if  9  is 
a  phase,  then  a  good  distortion  metric  might  be 

1 0(6-9)  =  1-cos  (9-0).  (3.1) 

Once  a  distortion  measure  is  chosen,  a  problem  of  interest  is  to  find  fundamental,  global  lower 
bounds  on  the  minimum  distortion  achievable  by  an  estimator  for  9.  Such  bounds  serve  as  bench¬ 
marks  for  evaluating  algorithms  that  approximate  estimators.  If  an  algorithm  achieves  the  global 
lower  bound,  then  it  cannot  be  improved  upon  in  terms  of  distortion  performance. 

Unfortunately,  past  research  in  this  area  is  incapable  of  dealing  with  periodic  distortion  functions 
such  as  (3.1).  However,  while  pursuing  our  primary  research  goals,  we  developed  the  following  [51]: 

Result  8  A  global  lower  bound  on  estimation  error  for  vector  parameters  where  some  of  the  dis¬ 
tortion  functions  are  periodic. 

s  • 

B.  Stability  of  Nonlinear  Least  Squares  Problems 

Another  result  we  have  established  while  pursuing  our  primary  goals  is  in  regards  to  the  stability 
of  nonlinear  least  squares  problems.  A  fair  number  of  problems  in  signal  and  image  processing  can 
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be  reduced  to  the  solution  of  a  nonlinear  least  squares  problem: 

y  =  argmin  ||/(y)  -  b||2  (3.2) 

y 

where  y  6  Rm,  b  6  Rn,  and  /  :  Rm  »->  Rn  is  a  nonlinear,  smooth  (i.e.  twice  continuously 
differentiable)  function. 

Often,  when  considering  stability  of  (3.2),  one  is  forced  to  consider  the  problem  in  a  stochastic 
setting.  There,  one  assumes  that  the  measurements  are  Gaussian  distributed  with  known  variance, 
b  ~  N(f(y),  cr2/),  where  I  is  the  n  x  n  identity  matrix.  With  this  assumption,  it  then  becomes 
possible  to  compute  Cramer-Rao  bounds  on  variance  of  unbiased  estimators  of  y  given  6.  However, 
the  resulting  bounds  are  valid  only  if  the  measurements  are  Gaussian.  In  the  presence  of  small, 
systematic  errors,  the  measurement  errors  may  be  better  modeled  as  deterministic,  instead  of 
stochastic.  In  [49],  we  proved  the  following  result: 

Result  9  There  are  a  set  of  conditions  which  guarantee  stability  of  (3.2)  in  both  the  traditional 
stochastic  sense,  and  in  the  deterministic  sense.  Hence ,  if  these  conditions  are  satisfied ,  then  the 
Cramer-Rao  lower  bound  on  the  variance  of  unbiased  estimators  is  finite,  and  the  estimate  y  varies 
smoothly  with  any  (deterministic)  perturbation  of  the  data  b ,  provided  the  perturbation  is  sufficiently 
small. 

We  have  applied  this  result  to  two  classic  problems  in  signal  processing  with  the  following  new 
results: 

Result  10  The  sinusoid  retrieval  problem  is  deterministically  stable. 

Result  11  The  multichannel,  blind  deconvolution  problem  is  deterministically  stable. 
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Chapter  4 

Conclusions 


In  this  report,  we  have  outlined  the  research  we  have  performed  while  sponsored  by  the  JSEP 
fellowship.  We  have  made  significant  and  fundamental  contributions  to  the  areas  of  tomography  in 
the  presence  of  motion,  fast  algorithms  for  tomographic  reconstruction,  and  stability  and  estimation 
problems  in  signal  processing. 
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